function [NumPores, AveDiameter, fibril_content, AllDiameters] = PoreSize(image1, pixel_size, strel_disk, erode_cutoff, top_hat, weiner, filter_size)
%% COLLAGE PORE SIZE ANALYSIS
% Matthew Zanotelli
%Updated by Adam Munoz- Feb 2018 to work with R2017A
%Updated by PVT - July 2018

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Image load sequence
% Input calibration information
calibration = (1/pixel_size);

imagedouble = image1;

% Weiner adaptive noise filtering
if weiner == 1
    image2 = wiener2(image1,[filter_size filter_size]); 
    image3 = image2 - min(image2(:));
else
    image3 = image1;
end

% Top hat filtering
if top_hat == 1
    se = strel('disk',strel_disk);
    tophatFiltered = imtophat(image3,se);
    J = imadjust(tophatFiltered);
else
    J = image3;
end


upper = mean(imagedouble(:)) + 3*std(imagedouble(:));
lower = mean(imagedouble(:)) - 3*std(imagedouble(:));
if imagedouble > upper | imagedouble < lower;
    T=2;
else
    T=graythresh(J);
end

% Make image binary and apply median filter to remove remaining noise
image4 = max(J(:))-J; % inverts pixel intensities
image5 = image4-min(image4(:)); % subtract min pixel intensity from image
image5 = image5-mean(image5(:)); % subtract mean pixel intensity from image
BW = im2bw(image5,0.01); % generate BW image for pore size autocorrelation
image5 = BW;

BW3 = imcomplement(image5);
BW4 = bwareaopen(BW3, 5);
BW4 = imcomplement(BW4);
image5 = BW4;

fibril_content = 100*bwarea(BW3)/1024^2; % was BW3

% troubleshooting

Area = bwarea(image5);
TotEroded = 0;
for i = 1:15;
    Erode = imerode(image5,strel('disk',i));
    AreaErode = bwarea(Erode);
    Eroded = AreaErode/Area;
    if Eroded <= erode_cutoff;
        break
    end
end
Erode = bwareaopen(Erode, 10);
fibril_content = 100*bwarea(Erode)/1024^2;
%% Output displays
% Display original image
figure;
imshow(image1);
% Display filtered image
figure
imshow(image5)
% Display eroded image; image to be analyized
figure;
imshow(Erode);

%% Output values
stats = regionprops(Erode,'EquivDiameter');
AllDiameters = [stats.EquivDiameter];
NumPores = size(AllDiameters,2);
AveDiameter = mean(AllDiameters)*calibration

end
